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Abstract 

We discuss some contradictions found in the literature concerning the problem of stability of 
coUisionless spherical stellar systems which are the simplest anisotropic generalization of the well- 
known polytrope models. Their distribution function F{E, L) is a product of power-low functions 
of the energy E and the angular momentum L, i.e. F oc L~''{—E)''. On the one hand, calculation 
of the growth rates in the framework of linear stability theory and N-body simulations show that 
these systems become stable when the parameter s characterizing the velocity anisotropy of the 
stellar distribution is lower than some finite threshold value, s < Scrit- On the other hand Palmer 
& Papaloizou (1987) showed that the instability remained up to the isotropic limit s — 0. 

Using our method of determining the eigenmodes for stellar systems, we show that the growth 
rates in weakly radially-anisotropic systems are indeed positive, but decrease exponentially as the 
parameter s approaches zero, i.e. 7 cx exp( — s*/s). In fact, for the systems with finite lifetime 
this means stability. 



Keywords: Galaxy: center, galaxies: kinematics and dynamics. 

1 Introduction 

Stability properties of stellar spherical clusters determine a set of dynamically allowable equilibrium 
configurations. Presence of an instability can, for example, lead to ellipsoidal deformation. 

For a long time, one believed that any stellar spherical clusters, except for the pathological models, 
were stable. This belief appeared after the classical works by Antonov (1960, 1962) devoted to isotropic 
systems, and was reassured in the following papers concerning some particular anisotropic systems 
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(e.g., Mikhailovsky et al. 1970, Doremus et al. 1971). So, it needed some time to realize the very 
possibility of instability of spheres, starting from Polyachenko & Shukhman (1972) until Merritt & 
Aguilar (1985), Barnes et al. (1986) and May & Binney (1986) have made it widely known. 

The radial orbit instability is suppressed for sufficiently rounded orbits, or when the kinetic energy 
stored in transverse directions T± becomes sufficiently high. Polyachenko & Shukhman (1981) pro- 
posed a global anisotropy parameter as a ratio of the radial to transverse kinetic energy of the system, 
^ = 2Tr/Tj_. For the Idlis model (Idlis 1956), they found the stability boundary ^ = 1.59. Later, 
Fridman & Polyachenko (1984), using this and two other families of models, proposed a hypothesis 
that the global anisotropy parameter can give a general stability criterion for anisotropic systems; the 
system is stable if < 1.7 ± 0.250 

Following studies of spherical models by means of the linear stability analysis (Saha 1991, Weinberg 
1991, Bertin et al. 1994) and N-body simulations (Merritt & Aguilar 1985, Barnes et al. 1986, 
Merritt 1987, Dejonghe & Merritt 1988, Meza & Zamorano 1997) included a variety of models radially 
anisotropic on the periphery and isotropic in the center, and vice versa. In these works, the stability 
boundaries in terms of the global anisotropy parameter fall in the broad range 1.2 < ^ < 2.9. Thus, 
the hypothesis about universal stabilization in the narrow region of ^ was denied. Note, however, that 
in each case a certain value of the stability boundary corresponding to the radially anisotropic system 
(^ > 1) was found. The same is referred to the generalized polytropes, with the DF 

F{E,L)^C{s,q)L-%-Ey , (1.1) 

which become stable at ^ sa 1.4 (Fridman & Polyachenko 1984; Barnes et al. 1986). Here E and L 
denote the energy and the angular momentum of a star, C(s, q) is the normalizing constant, s and q 
- parameters of the model. Additive constant in gravitational potential $0 (f) is chosen in such a way 
that $o(-R) = 0, where R is the radius of the system. 

Generalized polytropes are the simplest generalization of the isotropic polytrope models (the latter 
correspond to s = 0). The polytrope models are classical ones in stability theory of both gaseous 
and coUisionless gravitating systems. One can recall, for example, the work by Antonov (1962), in 
which stability of polytrope models with decreasing DF was shown. Models with increasing DF can 
be unstable, but they give exotic mass distribution with increasing density outwards, thus describing 
unrealistic stellar systems. 

The generalized polytropes are more versatile. The global anisotropy parameter for (jl.ip can be 
obtained in a simple form: 



Note that for this model the (local) anisotropy parameter /3(r) (Binney 1980) does not depend on 
radius, /S — 1 — 1/^ = s/2. It is possible to evaluate expressions for the radial and the transverse 
kinetic energy when s < 2, the limit s — > 2 corresponding to the system in which almost all orbits 
are radial. Since s = case is stable (let us consider only realistic models with q > 0), and the case 
s — >■ 2 is unstable due to the radial orbit instability, there should be a critical value of the parameter 
s — Scrit which divides stable and unstable systems. 

^Results of stability analysis of one of three families of models described in Fridman & Polyachenko (1984) were 
reconsidered later in Polyachenko (1987) report. Reconsidered stability boundary 2Tr/T±^ fell between 2.05 and 2.10, 
instead of the previous boundary, 2Tr/Tj_ = 1.62, i. e. the systems proved to be more stable than it was supposed 
before. The same boundary for this family was obtained later by Dejonghe and Merritt (1988) with the help of N-body 
simulations. A new (corrected) boundary is slightly out of the range 2Tr/Tj_ = 1.7 ± 0.25 suggested by Fridman & 
Polyachenko (1984). 
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Using the matrix method for spheres (Polyachenko & Shukhman 1981), which is analogous to 
the Kalnajs matrix method for disks (Kalnajs 1977), it was found that growth rates of the instabihty 
became small for s < 0.6, almost independently of parameter q (Fridman & Polyachenko 1984). Thus, 
the critical parameter for generalized polytropes is Scrit ~ 0.6 (or ^ w 1.4). Similar result was obtained 
by N-body simulations (Barnes et al. 1986). 

Palmer & Papaloizou (1987) (henceforth PP87) have investigated the same models using approx- 
imate equation for unstable modes with low growth rates. They showed that instability must persist 
even for models arbitrary close to isotropic limit s = 0; this seemingly contradicts previous results 
mentioned above. 

Solutions of the approximate equation form a set of infinite number of unstable modes with de- 
creasing growth rates (eigenvalues) which accumulate near zero frequency uj = 0. These eigenvalues 
correspond to eigenfunctions with different number of nodes; the nodeless eigenfunction gives the 
largest growth rate. As was argued in PP87, the largest eigenvalue cannot be caught by the approxi- 
mate equation, since it is too large to fulfill the assumed condition. 

Guided by their models B and C, which correspond to s = 2/3 and s = 1/3 {q = 1), we have 
calculated the growth rates near the isotropic limit using the approximate equation. The result was 
paradoxical: the largest eigenvalue kept on grow to infinity as s — )• 0, while one expects that all modes, 
including the nodeless one, would cease to zero. 

This paper pursues two goals. First, we try to reconcile the results obtained by the matrix method 
and N-body from one side, and the results of PP87 from the other side. Second, we clarify the paradox 
about the applicability of the approximate equation by PP87. For simplicity, we shall assume models 
with s < 1 only. This condition provides sufficiently smooth gravitational potential in the center, and 
linear dependence of the precession velocity on angular momentum for nearly radial orbits (for more 
details, see below). In Section 2, we derive an approximate equation for modes with low growth rates 
from the full integral equations for spheres obtained by Polyachenko et al. (2007). It coincides in all 
but one important detail with the approximate integro-differential equation by PP87 (the equivalence 
of two equations is demonstrated in Appendix). Numerical results are given in Section 3, Section 4 
contains conclusions. 

2 The integral equation for modes with low growth rates 

Traditional linear stability theories employ matrix methods, expanding the perturbed potential and 

density in scries using special biorthonormal sets of basis functions (Kalnajs 1977, Polyachenko & 
Shukhman 1981). As a result, one obtains a set of integral equations which incorporates the mode 
frequency in a complicated nonlinear manner. Thus each frequency is to be obtained separately by, 
for example, the Cauchy integration in the complex plane. 

Recently we have proposed an alternative method for calculation of eigenmodes (Polyachenko 2004, 
2005; Polyachenko et al. 2007). The advantages of our method are (i) linear form of the equation for 
eigenmodes and (ii) absence of the basic biorthonormal set that should be customized for a particular 
problem. The alternative method is the most adequate to derive the approximate integral equation 
similar to one used in PP87. We start with the full integral equation for perturbations proportional 
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to spherical harmonic with the index I: 



AttG ^ ^ r^i', f f dE' L' dV 



I 



X 



X U, i,(E,L;E',L') ' ' /'^ \ ' . (2.1) 

Integration in p.ip is over the curved triangle in the phase plane {E' , L'): $o(0) < i?' < 0, < L' < 
LciTc{E')] Lcirc(-E') is the angular momentum on the circular orbit with energy E; w is the eigenfre- 
quency; i}i^i^{E,L) = lini{E,L) + l2^2{E,L)] Vli^2 are the orbital frequencies; Vii^j'^F{E' , L') = 
Qi'^i'^{E', L') {dF /dE') + {dF/dL'); the coefficients are equal to zero for odd |/ — fc|, otherwise 

^ {i + ky.{i-ky. 

li and I2 are indices of expansion over angular variables wi and W2 in the action -angle formalism 
(Landau & Lifshitz 1976), 

S'^{h,l2,l3,Wi,W2) = ^{S<i>)i,i^{l) exp[j(;iwi + I2W2)], 
conjugated to the action variables li {i = 1, 2, 3): 



/i = i j \j2E- 2$o(r) - ^ dr, h = L, L,. 

Due to degeneracy on the azimuthal number m, the eigenfrequency w can be calculated for axially 
symmetric perturbations 6^{r,d;t) = x{'f)Pi{0)e~'^'^* only. The kernel of the integral equation is 

Tli,A,;i[A', {E, L- E', L') ^jdw^j dw[ Ti [r{E, L; w,), r'{E', L'; w[)] x 

X cosei,ijE,L,wi) cosei>j,^{E',L',w[), (2.2) 
where Qi^i^iE, L;wi) = {fti^i^/ili) wi - hScpiE, L;wi), 

r{E,L,wi) 

Sip{E, L,Wi) = L 



dx 



Ti(r,r') — r^^/r^-^^, ?'< = min(r, r'), r> = max(r, r'). 

Finally, the eigenfunctions 4>i^i^{E,L) are connected to the radial part of the perturbed potential as 
follows: 

TT 

4>iii2{E,L) = ^ J cosQij.^{E,L;wi)x[r{E,L,wi)] dwi. 
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To obtain the approximate equation by PP87 from (I2.ip (rather its full equivalent in the form of 
the integral equation in i5-space), one should make two simplifications, considering (i) low frequencies 
Lu = i"f and even spherical numbers I; {ii) domination of nearly radial orbits. 

The denominators of resonance terms l[ — I2 contain construction proportional to the preces- 
sion rate 

Qy^i-^ {E', L') ^i'2{n2-\n^) = I'^npAE', l'), (2.3) 

which is small for nearly radial orbits H 

np,{E,L) = w{E)L. 

Dropping the nonresonance terms and denoting 4>-^i^ hi^^^) = ^{E), from ()2.ip one can have 

dF 

For DF in the form F{E,L) — g{E) L^'', the approximate equation reads 

fe-2 Q 

where i'{E) = fli{E,0), Il{E,E') is the result of reduction of Ili^j^.i'^j'_^{E, L; E' , L') for the radial 
orbits. 

Due to singularity of DF, the integral in (12. 5p diverge when 7 = 0, and is large when 7 is small. 
This justifies omission of nonresonance terms, and sets constrains on the maximum value of 7. 

It is clear that main contribution to the integral comes from a narrow region L ~ 7/cc7, thus one 
can change the variable of integration and replace the upper boundary by infinity: 

l-^kwy-^ I{s), (2.6) 



7^ + ^2^2 L2 



where 





By appropriate change of the eigenfunction, one can reduce the problem to the integral equation 



X{s)-9{E)^ J dE'ns{E,E')^{E'), (2.8) 

-1 



^For generalized polytropes, the gravitational potential behaves like 'I"o(»') oc r'^~'' near the center. Thus, for the 
case of our interest s < 1 the gravitational force — <I'q(t") is non-singular at the center and hence the precession velocity 
of nearly radial orbits is indeed linear with respect to the angular momentum, ilp,-{E,L) ■uj{E) L (see, e.g., Touma 
& Tremaine 1997). Note that for singular 3>Q(r), (say, 3>o(?') oc r^, with p < 1) the dependence of precession rate on L 
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with 



A = 7^ (2.9) 
and the positively defined symmetric kernel function 



TZs{E,E')^a{s)^h{E)h{E') Q{E,E'), (2.10) 
where h{E) = gs{E)m''-'^{E)v{E), 

QiE.E')= j j , "-"-'■"■'7"'' (2.1, 

l2E' - 2$o(r') \ 2E - 2$o(?') 



^ (l^r)^ V D^k\ (2.12) 

^ ^ (2; + 1) sin(^s/2) ' ^ ^ 



In Appendix, we show that Eg. (|2.8p in i?-space is fully equivalent to the approximate integral 
equation in r-space obtained by PP87. 

The kernel (|2.10p defines a self-adjoint Hilbert-Schmidt operator in the infinite-dimensional space, 
so the eigenvalues A„ {n = 0,1,2...) must have an accumulation point, lim An = 0. Existence of 

n— >oo 

arbitrary small eigenvalues is crucial for PP87 in demonstrating the instability of singular generalized 
polytropes with s > 0. 

Note that in the limit s ^ 1 Eq. (|2.8I) with the kernel (I2.10p looks unnatural. Let us consider 
explicitly the case s = 0. The integral equation and the kernel then read as 



A-^{E) = J dE'no{E,E')^{E'), (2.13) 



■7? (P goiE)goiE')iy{E)iy{E') 

no{E,E)^a{0)^ ^^y^^ QiE,E), (2.14) 

with A = A(0), a(0) [327tG/{21 + 1)] J2 ■ For example, in the units where 4ttG = 1, for 

fe=2 

I = 2 one has a(0) = 3/5. A norm of the kernel is of order unity and thus first several eigenvalues, 
corresponding to eigenfunctions with few nodes, must be of order unity. It is needed to emphasize 
that there is no small parameter left in the problem (|2.13l) , the only small parameter s in the isotropic 
limit has disappeared from the equations. 

We have calculated several largest eigenvalues A„ for spherical indices I — 2,4,6. The results are 
summarized in the Table. The eigenvalues exceeding 1 are emphasized by boldface. In particular, for 
I ~ 2 two eigenvalues are greater than 1. This means that an arbitrary small anisotropy (or arbitrary 
small s) will produce exponentially high growth rates: 

7„ = A^" cx exp(i In A„) , n = 0, 1. (2.15) 



6 



l\n 





1 


2 


3 


4 


5 


2 


3.7851 


1.0301 


0.4762 


0.2637 


0.1621 


0.1068 


4 


1.3838 


0.4194 


0.2215 


0.1370 


0.0921 


0.0654 


6 


0.7029 


0.2209 


0.1232 


0.0803 


0.0566 


0.0419 



Table 1: The largest 6 eigenvalues A„ (i.e. A„ at s = 0) for quadrupole I = 2, and next two spherical harmonics {1 = 4 
and 1 = 5) for parameter q = 1.0. 



However, the growth rates of other modes are exponentiahy small: 

7„(5) = A„i/^(xexp(-i ln-1), n>2. (2.16) 

Note that (|2.15l) contradicts Eq. (12. 4p . in which the kernel becomes zero at s = due to the 
term dF/dL. The inconsistency evidently comes from changing the upper limit of integration in (|2.6p 



to infinity: for s — this integral turns into /q^""*" LdL/ (7^ + k'^tu'^ L^) and diverges if imax — >■ 00. 
However, such a form of the integrand is valid for nearly radial orbits only. Besides, we have expanded 
the integration region up to infinite angular momentum. These are justified for the systems mainly 
populated by nearly radial orbits, but not for the nearly isotropic ones. 

To cope with anomalously growing modes (|2.15p . one can take into account the finite value of 
imax in (|2.6|) . Changing the variable of integration, L — [7/ (fctu)]^, the integral (|2.6p can be reduced 
to 

L-'+^dL _ _2 f x~'+^ dx 



Since wLdrc ~ ^2{E,Lch-c{E)) — ^ni{E,LciTciE)) — il^ ■^k, where and >c are the circular and 
radial frequencies, one can replace the upper boundary of integration by r2/7, with some characteristic 
dynamical frequency fl ~ fin, flo = (G'A//(2i?^))^/^. Then instead of (|2.8p . we have 

n/7 ^_ 1 

/>,7/f^)= / ^^ = i / z^/^-i(l-zr/^dz (2.17) 
72/02 

(parameter j/fl ^ 1). At small s, this integral is finite, 

so /(0,7/f7) = ln(57/7), and the kernel turns to zero at s = 0. Unfortunately, it is impossible to 
take the integral explicitly for arbitrary s. Instead, in our crude approximation, we shall use a model 
expression obeying necessary features: 



%"/^)=^7— 7T— T =:T— 7T— T l-cxp{-.lnl^) . (2.18) 



2sin(2 7rs) L Vil/ J 2sin(2 7rs) 
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For very small s, s ^ [ln(r2/7)] ^ ^ 1, the expression gives /(s,7/ri) ~ ln(il/7), for s » [ln(r2/7) 
it coincides with old expression (|2.7p . 
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The corrected approximate integral equation can be written in the same form as the old one, 
(|2.10l) . but instead of (I2.9p one should write a new relation between A and 7: 

A^ . 7\^,s , 7=( , , tr.. Y- (2-19) 



i-(7/r!)" ' \i + x/n 

This relation clearly provides exponentially small growth rates for all modes at s — > 0. In this limit, 
the spectrum of eigenvalues A„ (s) remains the same, but now 



A, 



l/s 



asymptotically tends to zero when s — ^ 0. 

3 Unstable modes of generalized polytropes 

For numerical evaluations of growth rates of anisotropic polytropes, one usually introduces the fol- 
lowing units: 

47rG=l, A(s,(j) = 1, *(0) = 1, (3.1) 

where G is the gravitational constant, ^(r) is the relative potential, ^'(r) = — $o(r), and A{s,q) is 
the coefficient in the expression for density distribution, 

Po{r)^A{s,q)r-'^''+^. (3.2) 
The relation between A{s,q) and C{s,q) is well-known (see, e.g., Fridman & Polyachenko 1984): 

r((7 + i)r(-is + i) 



A{s,q) = C{s,q) (27r)3/2 2^^^/ 



2 



r(g + i(5-.)) ' 



r(x) denotes the Gamma function. 

In these units, the dimensionless potential ^/'(r) = — $o(^) satisfies the Poisson equation: 

(f->p 2 dip ^ „, 3-. , ^ 

dr'' r dr 

with a boundary condition ■0(0) = 1- To find a second boundary condition, one can notice that the 
Poisson equation possesses a solution in terms of a power series in z = r^"" = (see also Henon 
1973). Denoting ^/'(z) — ip{r), one has 

az^ 2 — s z dz (2 — s)^ z 

In the limit z < 1, the solution is ■(/5(z) = 1 + Diz + 0{z'^), where Di = -l/[(2 - s)(3 - s)]. Thus, 
one can obtain the potential by integrating p.4p from z = 0, with the boundary conditions iplO) — 1, 
(d'i/'/dz) = Di. The right boundary of integration z/? is determined from the condition ip{zji) — 0, 

thus radius of the system is i? = ^r^^ 
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Evaluation of the radial frequency ^{E) = D,i{E,0) is usual. The precession rate at low angular 
momenta (n7(i?) — {dilpr{E , L) / dL)\L=o) is calculated according to Eq. (2.11) by Polyachenko et al. 
(2007): 



m{E) 



v{E) 



dr 



1 



1 



2E + 2^{r) J2E + 2 



{E)J2E + 2 



The following asymptotic formulas, applicable for stars in the very center of the sphere, e = E + 1 <C 1, 
are used for testing of the precision of calculations: 



• s < 1 



• otherwise, s < 1 



1 / 262 

i+ieoH^--'^ 



s 

12^ 



r(- + -) 111 
^' r(i + i) 

r(i + i)r(i-i) 2 2 
r(i + i)r(i-i) 



(3.5) 
(3.6) 



where = p{l+p),p = 2 — s. 

Some problems were experienced in calculating the function Q{E,E') ()2.1ip . which is a part of 
the kernel (|2.10p . The following asymptotic formulas are useful (valid for I ~ 2): 



• on the diagonal e = £'<Cl,s = 

Q(e,e)c.y|(i + G)£-i/2 = 1.734 

where G = 0.915 965 594... is the Catalan's constant; 

• e' <C e <C 1, arbitrary s 



/2 



Q(.,.') = ^j|^i^-.-(e')*-. (3.7) 

The integral equation (j2.8p have been solved for modes with the spherical harmonics I = 2,4,6. 
First of all, we were interested in models with small parameter s ^ 1 to find out how small the growth 
rates of weakly anisotropic systems are. Also we have calculated growth rates for models B (s = 2/3, 
(7=1) and C (s = 1/3, q = 1) of PP87. In the numerical work, we restricted ourselves by models 
within the range < s < 0.8. 

Fig. [T] shows the behavior of the characteristic frequency ilois) = {GM/2R^y^^ for g = 1 and 
q — 0.7. One can see that the dependence is weak and fto ^ 0.1 in the considered range of the 
parameter s. 
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0.12 



0.1 - ~ ~ - - - - 




0.02 - 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

S 

Figure 1: The dependence of characteristic frequency 57/) on the parameter s for q ~ 1 (sohd hne) 
and q = 0.7 (dashed hne). 




Figure 2: The dependence of the scaled growth rates (7„ = 7„/J7d of four most unstable modes on 
the parameter s for the model q = I: dashed curves show solutions for (12.91) . solid curves - for the 
corrected relation (12.191) . Crosses show values of tin for model B, stars - for model C from PP87. The 
characteristic frequency is 57 = 0.08. 

Fig. [5] shows the growth rates <t„ in units of characteristic frequency flo v.s. the parameter s for 
q ~ 1. Dashed lines show growth rates obtained from (j2.8p with relation (12. 9^ : this case is equivalent 
to the approximate equation used by PP87. For s = 1/3 and s = 2/3, our growth rates agree 
satisfactorily with one obtained in the cited paper. Numbers denote modes (0 - the nodeless mode, 
1 - the mode with one node, etc.). The growth rates of the first two modes increase violently as the 
model approach the isotropic limit, other modes decrease exponentially to zero. 

Solid lines show the growth rates v.s. the parameter s obtained from (|2.19p . Its behavior complies 
with intuitive expectations that the unstable modes should be stabilized in the isotropic limit. 

A region near s « is magnified in Fig.[3l Due to fast (exponential) decrease of the growth rates, 
they become negligibly small at s = 0.01. 

Fig.lHshows several first eigenfunctions 5'„(£') {n — 0, 1, 2, 3) of the integral equation ()2.8p for the 
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Figure 3: The same as in Fig.[2l in Ig - Ig axis. 



model s = 1/3, q = I and corresponding radial parts of the potential x„(r): 

r^^AE) 

Xn{r)= h^hm^^niE) f f''^'^'^''^ (3.8) 

(for derivation of this relation, see Appendix). In all cases, the eigenfunctions have equal number of 
nodes coinciding with n. The form of radial parts of the potential Xn{f) is in qualitative agreement 
with those presented in Fig. 1 of PP87, which are the solutions of the integral equation (2.25). 

It is interesting to compare the growth rates obtained with our approximate integral equation with 
independent calculations of generalized polytropes. In Fig.[5l we show the dependence of growth rates 
for first seven modes on the parameter s for q = 0.7. Crosses mark an "experimental" curve from 
Fridman & Polyachenko (1984). The curve breaks at s « 0.6, 7 « 0.004 because of the accuracy of 
the matrix method employed. The numbers poorly agree with each other, so that it is impossible to 
join the curve marked by crosses with a new curve of the principal mode n — Q. 

The validity of the approximate integral equation is restricted to very small growth rates. For 
modes with sufficiently high n, this restriction does not play any role since their growth rates are small 
for all s, but it is important for modes with few nodes, especially for the nodeless one. Its growth 
rate increases with anisotropy, and presumably at s < 0.15 achieves the upper boundary. This might 
be the reason of the evident discrepancy between the two curves. In Fig.[5l the uncertainty region 
0.15 < s < 0.6 is shown by a dashed line. 



4 Conclusions 

1. It is difficult to reconcile the results of the linear stability analysis (Fridman & Polyachenko, 1984) 
and N-body experiments (Barnes et al. 1986) with the results by PP87 for two reasons. 

First, the approximate integral equation derived in PP87, on which their analysis is based, is 
applicable to very small growth rates only. In the isotropic limit corresponding s — ^ 0, these growth 
rates are exponentially small, i.e. 7 cx exp(— s,/s). The estimates show, that even for s ~ 0.5, the 
allowed growth rates are much less than O.OlfJ^), that is definitely below any reasonable accuracy of 
the matrix method and N-body experiments. 
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1 2 3 4 5 6 

r 

Figure 4: The first four eigenfunctions of quadrupole harmonics I — 2: 5'„(-E) (upper figure) and the 
radial part of the potential Xn{r) (lower figure). The normalization of the eigenfunction is arbitrary. 




(1984). The characteristic frequency is taken = 0.02. 
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Second, strictly speaking, the approximate integral equation is not applicable to unstable modes 
having eigenfunctions with just few nodes, including the principal mode with maximum growth rate, 
since for not too large n the eigenvalues A„ are not very small even in the limit s — 0. However the 
principal mode plays the major role in determining the stability boundary (see Fridman & Polyachenko 
(1984), Barnes et al. (1986)). 

2. The growth rates formally accurate when they are small. However, it is likely that once 
an = jn/^D < 1, they give reasonable estimates of actual growth rates in practice (PP87). Then 
since isotropic systems with decreasing DF are stable (Antonov 1960, 1962), all modes should become 
stable when s approaching zero, and valid approximate equation must describe all modes correctly. 
This fact is in contradiction to our solution of PP87's approximate equation (2.25) (equivalent of our 
Eq. (|2.8p in which the relation A = 7** is assumed). This solution demonstrates explicitly that for the 
quadrupole harmonic I — 2 there are two modes with exponentially increasing growth rates for s — ?> 0. 

The reason for such a discrepancy arises from behavior of the terms of the approximate integral 
equation in the limit s ^ 0. Using the notation of PP87, the growth rates can be expressed in the 
form 

^ A1-A2' 

where Ai are some averages (quadratic forms) of the positively defined operators]! In particular, A3 
denotes the average of the integral operator defined by p.lOp . A2 denotes the omitted nonresonance 
part and Ai denotes the operator which is a left side of the radial Poisson equation. Since both terms 
Al and A3 retain in the limit when the system becomes isotropic. Palmer & Papaloizou infer that 
instability exists no matter how weak the divergence in DF as L — >■ is. 

In this paper we argue that the term A3 must vanish when s ^ in order to comply with sta- 
bilization of isotropic models. Using our alternative method of determining the unstable eigenmodes 
based on solution of the linear eigenvalue problem, we derived the appropriate integral equation. This 
equation gives a set of unstable modes, all of which become stable in the isotropic limit s — > 0. 

To summarize, our considerations prove that the instability growth rates of all modes in unbounded 
models at L = indeed do not vanish unless the models are isotropic, but they becomes exponentially 
small. Actually it means stability if we take into account a finite lifetime of real astronomical objects. 
Besides, the most probable distributions are non-singular ones, so we have to infer that stable dis- 
tributions generally become unstable at some finite value of radial anisotropy, i.e. finite anisotropy 
threshold exists. Width of the threshold depends on a particular model. 
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APPENDIX A. Equivalence of the integral equation in £'-space and the 
integro-difFerential equation (2.25) of PP87 in r-space 

Actually, the integro-differential equation (2.25) of PP87 presents the radial part of the Poisson 
equation for Z-th harmonic of the potential 



1 d ^2 dx{r) Ijl + 1) 
dr dr 



r^^-^x(r)=4.Gn(r), {Al) 



where n(r) is the radial part of a density perturbation Sp{r,9;t) = n(r)P;(cos0) e~'"*. Wc can 
calculate it using an expression for the perturbed DF Sf obtained from the kinetic equation written in 
the actions - angles variables (Landau & Lifshitz 1976). For details of calculation of 6f see Polyachenko 
& Shukhman (1981) or Fridman & Polyachenko (1984). Thus we have 

n(r)P;(cos6') = j dv5f{r,v) = 
= -EE / dvPl' (O)^r'^ (Bin go)e"--0M. (E, L) ^^'^^^ + .HW^.+l.^.) ^ (^2) 

11 (2 

where sin^o = L^/L. Separation of the radial part of perturbed density yields: 

n(r) = \{2l + l) j Sp{r,0)Pi{cose) smOdO. (A3) 



We can integrate over using addition theorem for Legendre polynomials and a relation which connects 
6 with the angular variable W2 

cos 6' = cos 00 cos{w2 — dSi/dh), {Ai) 

where Si is radial action, 



5i(r)= / dr' ^/2E - 2^o{r') - L'^ /r''^ . 



We have 

Pi{cose)= e-^'=("'=-^^i/^^^)P;'=(sin0o)Pr''(O)e"'''^- (^5) 



k=-l 

Since 



, OA ^ ^ /i sill ^0 L smOo 
Vr = ±\ 2E-2<^o{r)- V0=±-\1 = 

then we obtain for the volume element in the velocity space 

ALdLdEdisiaeo) 



dv = dvrdvedvu, = 



r2y^2£-2$o(r-)-LV^^ y^cos^ fi^o - cos^ 61 
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(The factor 4 appears here because we have to take into account particles having velocities Vr and vo 
of both signs.) To integrate over 9 in (A3) it is convenient to go to integration over W2- We have 

sin^d^ = cos^o sm{w2 — dSi/dh) dw2 

and 

271- 

{A6) 



^ J cos^ ^0 — cos^ 6 ^ 



As a result we obtain 



n(r) = -l(2^ + l)(2.)EE / / 4L dLdE f ' ^) ^^--'^''^'^'^ h OF/81. + I2 BFBI^ ^ 



1 

X 

-1 

Taking into account that the integral over z is 

1 



j dzP}^ (0) {z) Pf^^ (0) P/^ {z), z = sin 60- {A7) 



J dz P/^ (O)Pr'^ (^)Pr'^ (0)P/^ (2) = ^ Dl^ , 



we find the final expression for the radial part of perturbed density n(r): 



ijj ^2£-2$o(r)-LVr2 " ^^W. 



r 

(l (2 = 



where (see Sec. 2) 

111 /■ rfrx(?')coseii;2 



2£;-2$o(r)-iVr2 



Keeping in the sum over li the resonance sunimand li = —^h and supposing growth rate to be small, 
we retain the contribution with dF/ dL only: 

~ -^l^D^H^ LdL cose_j,/2,,,. (A9) 

Since the leading contribution to the integral over L comes from small L, we put where it is possible, 
L = and suppose Opr = va{E) L. Then @_i^/2,i2{E,0) = h^^, and finally 



^%t^2 ^2E-2Mr) [ ^ / 



r„ax(i5) 

v{E) f dr'x{r') 



2E-2^o{r') 
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Substitution to r.h.s. of (Al) yields integro-differential equation (2.25) of PP87 in r-space: 

R. 

|:r2^-U^ + l)xM = -^ j K{r,r')x{r')dr', [All) 



with the kernel 

sin(7rs/2)^ J JoK - 0<bJ^\ . loK - 0<bJ^'\ 



k=2 r .j2E-2^o{r) j2E-2<^o{r') 

max[$o(r),*o(i-')J * * 

Now it is easy to demonstrate that (All) is equivalent to the integral equation (2.8) in £'-space. We 
write (AlO) in the form 

47r2s 1 1 ^ .,^rg{E)uj^-\E)dE 



sm(7rs/2) 7^ r- J 2E - 2Mr) 



where 

^ { ^2E-2Mr') 

is the radial part of perturbed potential x(r) averaged over radial orbit. Rewriting the Poisson 
equation (Al) in the integral form, 

X{r) = j dr'r'^W)Mry), (A15) 

and averaging both part of this equation over radial orbit with the energy E according to (A14), we 
obtain the integral equation 

^(E) = — / dE'JC(E,E')^(E'), (A16) 
Y J 



]C{E,E') = ^^p^ g{E')w'-\E')v{E)x 



where 

rCK F,'\ = 

(2Z + 1) sin(7rs/2) 



X 



^Dffc" J dr J dr' 



k=2 i i ^2E'-2Mr') x/2E-2^o{r) 

Symmetrizing this equation with the help of the substitution 



we obtain the integral equation (2.8) with the kernel (2.10) presented in the main text. 

Finally, using the expression (A13) for n(r) and the integral form of the Poisson equation (A15), we 
can easily obtain the relation p.8p . which connects the eigenfunction 'i'n{E) of the integral equation 
(|2.8p with the eigenfunction x„(r) of PP87's integro-differential equation (2.25) (or Eq. (All) in our 
notations). 
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